cafe afternoon wait
1 1 0 4.9989926
2 1 1 2.2133944
3 1 0 4.1866730
4 1 1 3.5624399
5 1 0 3.9956779
6 1 1 2.8957176
7 1 0 3.7804582
8 1 1 2.3844837
9 1 0 3.8617982
10 1 1 2.5800004
11 2 0 2.7421223
12 2 1 1.3525907
13 2 0 2.5215095
14 2 1 0.9628102
15 2 0 1.9543977Introduction à la modélisation statistique bayésienne
Un cours en R, Stan, et brms
Ladislas Nalborczyk (LPC, LNC, CNRS, Aix-Marseille Univ)
Planning
Cours n°01 : Introduction à l’inférence bayésienne
Cours n°02 : Modèle Beta-Binomial
Cours n°03 : Introduction à brms, modèle de régression linéaire
Cours n°04 : Modèle de régression linéaire (suite)
Cours n°05 : Markov Chain Monte Carlo
Cours n°06 : Modèle linéaire généralisé
Cours n°07 : Comparaison de modèles
Cours n°08 : Modèles multi-niveaux
Cours n°09 : Modèles multi-niveaux généralisés
Cours n°10 : Data Hackathon
Modèles multi-niveaux
Le but est de construire un modèle qui puisse apprendre à plusieurs niveaux, qui puisse produire des estimations qui seront informées par les différents groupes présents dans les données. Nous allons suivre l’exemple suivant tout au long de ce cours.
Imaginons que nous ayons construit un robot visiteur de cafés, et que celui-ci s’amuse à mesurer le temps d’attente après avoir commandé. Ce robot visite 20 cafés différents, 5 fois le matin et 5 fois l’après-midi.
Robot et café
Robot et café
df %>%
ggplot(aes(x = factor(cafe), y = wait, fill = factor(afternoon) ) ) +
geom_dotplot(
stackdir = "center", binaxis = "y",
dotsize = 1, show.legend = FALSE
) +
geom_hline(yintercept = mean(df$wait), linetype = 3) +
facet_wrap(~afternoon, ncol = 2) +
labs(x = "Café", y = "Temps d'attente (en minutes)")
Robot et café, premier modèle
On peut construire un premier modèle, qui estime le temps moyen (sur tous les bistrots confondus) pour être servi.
\[ \begin{align} w_{i} &\sim \mathrm{Normal}(\mu_{i}, \sigma) \\ \mu_{i} &= \alpha \\ \alpha &\sim \mathrm{Normal}(5, 10) \\ \sigma &\sim \mathrm{HalfCauchy}(0, 2) \\ \end{align} \]
Half-Cauchy
\[ p(x\ |\ x_{0}, \gamma) = \Big(\pi \gamma \Big[1 + \big(\tfrac{x-x_{0}}{\gamma}\big)^{2}\Big] \Big)^{-1} \]
Robot et café, premier modèle
Diagnostic plot
Un intercept par café
Deuxième modèle qui estime un intercept par café. Équivalent à construire 20 dummy variables.
\[ \begin{align} w_{i} &\sim \mathrm{Normal}(\mu_{i}, \sigma) \\ \mu_{i} &= \alpha_{\text{café}[i]} \\ \alpha_{\text{café}[i]} &\sim \mathrm{Normal}(5, 10) \\ \sigma &\sim \mathrm{HalfCauchy}(0, 2) \\ \end{align} \]
Un intercept par café
Estimate Est.Error Q2.5 Q97.5
b_factorcafe1 3.4507647 0.25018675 2.9664917 3.9386539
b_factorcafe2 1.7307201 0.25374784 1.2485055 2.2264165
b_factorcafe3 3.3249216 0.26433504 2.8021550 3.8301851
b_factorcafe4 2.7974581 0.26425751 2.2751067 3.3239308
b_factorcafe5 1.4623576 0.26243793 0.9384310 1.9813639
b_factorcafe6 3.6400246 0.26410037 3.1291509 4.1543382
b_factorcafe7 2.9411357 0.26073639 2.4383438 3.4565863
b_factorcafe8 3.1704641 0.26487715 2.6407843 3.7010120
b_factorcafe9 3.3324917 0.27025623 2.8140962 3.8654680
b_factorcafe10 3.0991653 0.26583897 2.5846845 3.6267451
b_factorcafe11 1.9201884 0.25509560 1.4049837 2.4281218
b_factorcafe12 3.4904030 0.26155442 2.9813862 4.0094279
b_factorcafe13 3.2200183 0.25894415 2.7217235 3.7323040
b_factorcafe14 2.6283465 0.26416436 2.1115814 3.1708637
b_factorcafe15 3.4819173 0.26243431 2.9717212 3.9893268
b_factorcafe16 3.0037880 0.25855590 2.4799206 3.5094136
b_factorcafe17 3.8779458 0.26108680 3.3724571 4.3947832
b_factorcafe18 5.5257080 0.26203818 5.0029700 6.0331673
b_factorcafe19 2.9703972 0.26637172 2.4491741 3.5120111
b_factorcafe20 3.3632874 0.25660563 2.8579161 3.8496222
sigma 0.8218442 0.04495196 0.7409397 0.9170816
lprior -66.1597161 0.02866500 -66.2218965 -66.1075367
lp__ -310.3586844 3.50109680 -317.9728987 -304.6931867Modèle multi-niveaux
Est-ce qu’on ne pourrait pas faire en sorte que le temps mesuré au café 1 informe la mesure réalisée au café 2, et au café 3 ? Ainsi que le temps moyen pour être servi ? Nous allons apprendre le prior à partir des données…
\[ \begin{align} \text{Niveau 1}: w_{i} &\sim \mathrm{Normal}(\mu_{i}, \sigma) \\ \mu_{i} &= \alpha_{\text{café}[i]} \\ \text{Niveau 2}: \alpha_{\text{café}} &\sim \mathrm{Normal}(\alpha,\sigma_{\text{café}}) \\ \alpha &\sim \mathrm{Normal}(5, 10) \\ \sigma_{\text{café}} &\sim \mathrm{HalfCauchy}(0, 2) \\ \sigma &\sim \mathrm{HalfCauchy}(0, 2) \\ \end{align} \]
Le prior de l’intercept pour chaque café (\(\alpha_{\text{café}}\)) est maintenant fonction de deux paramètres (\(\alpha\) et \(\sigma_{\text{café}}\)). \(\alpha\) et \(\sigma_{\text{café}}\) sont appelés des hyper-paramètres, ce sont des paramètres pour des paramètres, et leurs priors sont appelés des hyperpriors. Il y a deux niveaux dans le modèle…
Équivalences (encore)
\[ \begin{align} w_{i} &\sim \mathrm{Normal}(\mu_{i}, \sigma) \\ \mu_{i} &= \alpha_{\text{café}[i]} \\ \alpha_{\text{café}} &\sim \mathrm{Normal}(\alpha, \sigma_{\text{café}}) \\ \end{align} \]
NB : \(\alpha\) est ici défini dans le prior de \(\alpha_{\text{café}}\) mais on pourrait, de la même manière, le définir dans le modèle linéaire :
\[ \begin{align} w_{i} &\sim \mathrm{Normal}(\mu_{i}, \sigma) \\ \mu_{i} &= \alpha + \alpha_{\text{café}[i]} \\ \alpha_{\text{café}} &\sim \mathrm{Normal}(0, \sigma_{\text{café}}) \\ \alpha &\sim \mathrm{Normal}(5, 10) \\ \end{align} \]
On peut toujours “enlever” la moyenne d’une distribution gaussienne et la considérer comme une constante plus une gaussienne centrée sur zéro.
NB : quand \(\alpha\) est défini dans le modèle linéaire, les \(\alpha_{\text{café}}\) représentent des déviations de l’intercept moyen. Il faut donc ajouter \(\alpha\) et \(\alpha_{\text{café}}\) pour obtenir le temps d’attente moyen par café…
Équivalences (encore)
Modèle multi-niveaux
Ce modèle a 23 paramètres, l’intercept général \(\alpha\), la variabilité résiduelle \(\sigma\), la variabilité entre les cafés \(\sigma_{\text{café}}\), et un intercept par café.
Shrinkage
Shrinkage magic (Efron & Morris, 1977)
The James-Stein estimator is defined as \(z = \bar{y} + c(y - \bar{y})\), where \(\bar{y}\) is the grand average, \(y\) an individual estimation and \(c\) a constant, the “shrinking factor” (Efron & Morris, 1977).
Shrinkage magic (Efron & Morris, 1977)
The James-Stein estimator is determined both by the variability in the measure (e.g., its standard deviation, influencing the shrinking factor \(c\)) and by its distance to the average estimation (i.e., \(y - \bar{y}\)). In other words, extreme observations are trusted less. Shrinkage therefore acts as a safeguard against overfitting in multilevel models.
Pooling
Le shrinkage observé slide précédente est dû à des phénomènes de partage (pooling) de l’information entre les cafés. L’estimation de l’intercept pour chaque café informe l’estimation de l’intercept des autres cafés, ainsi que l’estimation de l’intercept général (i.e., la moyenne générale).
On distingue en général trois perspectives (ou stratégies) :
Complete pooling : on suppose que le temps d’attente est invariant, on estime un intercept commun (
mod1)No pooling : on suppose que les temps d’attente de chaque café sont uniques et indépendants: on estime un intercept par café, mais sans informer le niveau supérieur (
mod2)Partial pooling : on utilise un prior adaptatif, comme dans l’exemple précédent (
mod3)
La stratégie complete pooling en général underfit les données (faibles capacités de prédiction) tandis que le stratégie no pooling revient à overfitter les données (faibles capacités de prédiction ici aussi). La stratégie partial pooling (celle des modèles multi-niveaux) permet de balancer underfitting et overfitting.
Comparaison de modèles
On peut comparer ces trois modèles en utilisant le WAIC (discuté au Cours n°07).
# calcul du WAIC et ajout du WAIC à chaque modèle
mod1 <- add_criterion(mod1, "waic")
mod2 <- add_criterion(mod2, "waic")
mod3 <- add_criterion(mod3, "waic")
# comparaison des WAIC de chaque modèle
w <- loo_compare(mod1, mod2, mod3, criterion = "waic")
print(w, simplify = FALSE) elpd_diff se_diff elpd_waic se_elpd_waic p_waic se_p_waic waic se_waic
mod3 0.0 0.0 -253.9 8.3 18.3 1.5 507.9 16.6
mod2 -0.9 1.3 -254.8 8.4 19.8 1.6 509.7 16.8
mod1 -57.2 10.6 -311.1 10.5 2.0 0.3 622.2 21.1 On remarque que le modèle 3 a seulement 18 effective parameters (pWAIC), et moins de paramètres que le modèle 2, alors qu’il en a en réalité 2 de plus… posterior_summary(mod3)[3, 1] nous donne le sigma du prior adaptatif des \(\alpha_{\text{café}}\) (\(\sigma_{\text{café}} = 0.82\)). On remarque que ce sigma est très faible et correspond à assigner un prior très contraignant, ou régularisateur.
Comparaison de modèles
On compare le premier modèle (complete pooling model) et le troisième modèle (partial pooling model).
Estimate Est.Error Q2.5 Q97.5
b_Intercept 3.119875 0.07985599 2.960233 3.274415
sigma 1.143102 0.05749243 1.035762 1.263518 Estimate Est.Error Q2.5 Q97.5
b_Intercept 3.1192355 0.20255593 2.7183344 3.5238007
sigma 0.8222354 0.04359112 0.7414241 0.9121809Les deux modèles font la même prédiction (en moyenne) pour \(\alpha\), mais le modèle 3 est plus incertain de sa prédiction que le modèle 1…
L’estimation de \(\sigma\) du modèle 3 est plus petite que celle du modèle 1 car le modèle 3 décompose la variabilité non expliquée en deux sources : la variabilité du temps d’attente entre les cafés \(\sigma_{\text{café}}\) et la variabilité résiduelle \(\sigma\).
Robot et café
Imaginons que notre robot ne visite pas tous les cafés le même nombre de fois (comme dans le cas précédent) mais qu’il visite plus souvent les cafés proches de chez lui…
Shrinkage
On observe que les cafés qui sont souvent visités (à droite) subissent moins l’effet du shrinkage. Leur estimation est moins “tirée” vers la moyenne générale que les estimations des cafés les moins souvent visités (à gauche).
Aparté : effets fixes et effets aléatoires
Cinq définitions (contradictoires) relevées par Gelman (2005).
- Fixed effects are constant across individuals, and random effects vary.
- Effects are fixed if they are interesting in themselves or random if there is interest in the underlying population.
- When a sample exhausts the population, the corresponding variable is fixed; when the sample is a small (i.e., negligible) part of the population the corresponding variable is random.
- If an effect is assumed to be a realized value of a random variable, it is called a random effect.
- Fixed effects are estimated using least squares (or, more generally, maximum likelihood) and random effects are estimated with shrinkage.
Gelman & Hill (2006) suggèrent plutôt l’utilisation des termes de constant effects et varying effects, et de toujours utiliser la modèlisation multi-niveaux, en considérant que ce qu’on appelle effet fixe peut simplement être considéré comme un effet aléatoire dont la variance serait égale à \(0\).
Régularisation et terminologie
Le fait de faire varier les intercepts de chaque café est simplement une autre manière de régulariser (de manière adaptative), c’est à dire de diminuer le poids accordé aux données dans l’estimation. Le modèle devient à même d’estimer à quel point les groupes (ici les cafés) sont différents, tout en estimant les caractéristiques de chaque café…
Différence entre les cross-classified (ou crossed) multilevel models et nested or hierarchical multilevel models. Le premier type de modèle concerne des données structurées selon deux (ou plus) facteurs aléatoires non nichés. Le deuxième type de modèles concerne des données structurées de manière hiérarchique (e.g., un élève dans une classe dans une école dans une ville…). Voir ce thread pour plus de détails.
Les deux types de modèles s’écrivent cependant de manière similaire, sur plusieurs niveaux. Le terme “multi-niveaux” (dans notre terminologie) fait donc référence à la structure du modèle, à sa spécification. À distinguer de la structure des données.
Exemple de modèle “cross-classified”
On pourrait se poser la question de savoir si la récence des cafés (leur âge) ne serait pas une source de variabilité non contrôlée ? Il suffit d’ajouter un intercept qui varie par âge, et de lui attribuer un prior adaptatif.
\[ \begin{align} w_{i} &\sim \mathrm{Normal}(\mu_{i}, \sigma) \\ \mu_{i} &= \alpha + \alpha_{\text{café}[i]} + \alpha_{\text{âge}[i]} \\ \alpha_{\text{café}} &\sim \mathrm{Normal}(5, \sigma_{\text{café}}) \\ \alpha_{\text{âge}} &\sim \mathrm{Normal}(5, \sigma_{\text{âge}}) \\ \alpha &\sim \mathrm{Normal}(0, 10) \\ \sigma_{\text{café}} &\sim \mathrm{HalfCauchy}(0, 2) \\ \sigma_{\text{âge}} &\sim \mathrm{HalfCauchy}(0, 2) \\ \sigma &\sim \mathrm{HalfCauchy}(0, 2) \\ \end{align} \]
Robot et café : varying intercept + varying slope
On s’intéresse maintenant à l’effet du moment de la journée sur le temps d’attente. Attend-on plus le matin, ou l’après-midi ?
\[ \begin{align} w_{i} &\sim \mathrm{Normal}(\mu_{i},\sigma) \\ \mu_{i} &= \alpha_{\text{café}[i]} + \beta_{\text{café}[i]} A_{i} \\ \end{align} \]
Où \(A_{i}\) est une dummy variable codée 0/1 pour le matin et l’après-midi et \(\beta_{\text{café}}\) est un paramètre de différence entre le matin et l’après-midi.
Remarque : on sait que les cafés ont des intercepts et des pentes qui covarient… Les cafés populaires seront surchargés le matin et beaucoup moins l’après-midi, résultant en une pente importante. Ces cafés auront aussi un temps d’attente moyen plus long (i.e., un intercept plus grand). Dans ces cafés, \(\alpha\) est grand et \(\beta\) est loin de zéro. À l’inverse, dans un café peu populaire, le temps d’attente sera faible, ainsi que la différence entre matin et après-midi.
On pourrait donc utiliser la co-variation entre intercept et pente pour faire de meilleures inférences. Autrement dit, faire en sorte que l’estimation de l’intercept informe celle de la pente, et réciproquement.
Robot et café : varying intercept + varying slope
On s’intéresse maintenant à l’effet du moment de la journée sur le temps d’attente. Attend-on plus le matin, ou l’après-midi ?
\[ \begin{align} w_{i} &\sim \mathrm{Normal}(\mu_{i},\sigma) \\ \mu_{i} &= \alpha_{\text{café}[i]} + \beta_{\text{café}[i]} A_{i} \\ \begin{bmatrix} \alpha_{\text{café}} \\ \beta_{\text{café}} \\ \end{bmatrix} &\sim \mathrm{MVNormal}\bigg(\begin{bmatrix} \alpha \\ \beta \end{bmatrix}, \textbf{S}\bigg) \\ \end{align} \]
La troisième ligne postule que chaque café a un intercept \(\alpha_{\text{café}}\) et une pente \(\beta_{\text{café}}\), définis par un prior Gaussien bivarié (i.e., à deux dimensions) ayant comme moyennes \(\alpha\) et \(\beta\) et comme matrice de covariance \(\textbf{S}\).
Aparté : distribution gaussienne multivariée
\[\mathbf{x} \sim \mathcal{N}(\boldsymbol{\mu}, \boldsymbol{\Sigma})\]
Où \(\boldsymbol{\mu}\) est un vecteur (à \(k\) dimensions) de moyennes, par exemple: mu <- c(a, b).
\(\boldsymbol{\Sigma}\) est une matrice de covariance de \(k \times k\) dimensions, et qui correspond à la matrice donnée par la fonction vcov().
\[ \begin{align} \boldsymbol{\Sigma} &= \begin{pmatrix} \sigma_{\alpha}^{2} & \sigma_{\alpha} \sigma_{\beta} \rho \\ \sigma_{\alpha} \sigma_{\beta} \rho & \sigma_{\beta}^{2} \\ \end{pmatrix} \\ \end{align} \]
Aparté : distribution gaussienne multivariée
Aparté : distribution gaussienne multivariée, calcul matriciel
\[ \begin{align} \boldsymbol{\Sigma} &= \begin{pmatrix} \sigma_{\alpha}^{2} & \sigma_{\alpha} \sigma_{\beta} \rho \\ \sigma_{\alpha} \sigma_{\beta} \rho & \sigma_{\beta}^{2} \\ \end{pmatrix} \\ \end{align} \]
Cette matrice peut se construire de deux manières différentes, strictement équivalentes.
Aparté : distribution gaussienne multivariée, calcul matriciel
\[ \begin{align} \boldsymbol{\Sigma} &= \begin{pmatrix} \sigma_{\alpha}^{2} & \sigma_{\alpha} \sigma_{\beta} \rho \\ \sigma_{\alpha} \sigma_{\beta} \rho & \sigma_{\beta}^{2} \\ \end{pmatrix} \\ \end{align} \]
La deuxième méthode est pratique car elle considère séparément les écart-types et les corrélations.
Robot et café : varying intercept + varying slope
\[ \begin{align} w_{i} &\sim \mathrm{Normal}(\mu_{i},\sigma) \\ \mu_{i} &= \alpha_{\text{café}[i]} + \beta_{\text{café}[i]} A_{i} \\ \begin{bmatrix} \alpha_{\text{café}} \\ \beta_{\text{café}} \\ \end{bmatrix} &\sim \mathrm{MVNormal}\bigg(\begin{bmatrix} \alpha \\ \beta \end{bmatrix}, \textbf{S}\bigg) \\ \textbf{S} &= \begin{pmatrix} \sigma_{\alpha} & 0 \\ 0 & \sigma_{\beta} \\ \end{pmatrix} \ \textbf{R} \begin{pmatrix} \sigma_{\alpha} & 0 \\ 0 & \sigma_{\beta} \\ \end{pmatrix} \\ \alpha &\sim \mathrm{Normal} (0, 10) \\ \beta &\sim \mathrm{Normal} (0, 10) \\ \sigma &\sim \mathrm{HalfCauchy} (0, 1) \\ \sigma_{\alpha} &\sim \mathrm{HalfCauchy} (0, 1) \\ \sigma_{\beta} &\sim \mathrm{HalfCauchy} (0, 1) \\ \textbf{R} &\sim \mathrm{LKJ}(2) \\ \end{align} \]
\(\textbf{S}\) est définie en factorisant \(\sigma_{\alpha}\), \(\sigma_{\beta}\), et la matrice de corrélation \(\textbf{R}\). La suite du modèle définit simplement les priors pour les effets constants. La dernière ligne spécifie le prior pour \(\textbf{R}\).
LKJ prior
D’après Lewandowski et al. (2009). Un seul paramètre \(\zeta\) spécifie la concentration de la distribution du coefficient de corrélation. Le prior \(\mathrm{LKJ}(2)\) définit un prior peu informatif pour \(\rho\) qui est sceptique des corrélations extrêmes (i.e., proches de \(-1\) ou \(1\)).
Rappels de syntaxe
Le package brms utilise la même syntaxe que les fonctions de base R (comme lm) ou que le package lme4.
La partie gauche représente notre variable dépendante (ou outcome, i.e., ce qu’on essaye de prédire).
Rappels de syntaxe
La première partie de la partie droite de la formule représente les effets constants (effets fixes), tandis que la seconde partie (entre parenthèses) représente les effets variables (effets aléatoires).
Le premier modèle ci-dessus contient seulement un intercept variable, qui varie par Subject. Le deuxième modèle contient également un intercept variable, mais aussi une pente variable pour l’effet de Days.
Rappels de syntaxe
Lorsqu’on inclut plusieurs effets variables (e.g., un intercept et une pente variables), brms postule qu’on souhaite aussi estimer la corrélation entre ces deux effets. Dans le cas contraire, on peut supprimer cette corrélation (i.e., la fixer à 0) en utilisant ||.
Modèle brms
On spécifie un intercept et une pente (pour l’effet d’afternoon) qui varient par cafe.
mod5 <- brm(
wait ~ 1 + afternoon + (1 + afternoon | cafe),
prior = c(
set_prior("normal(0, 10)", class = "Intercept"),
set_prior("normal(0, 10)", class = "b"),
set_prior("cauchy(0, 2)", class = "sigma"),
set_prior("cauchy(0, 2)", class = "sd")
),
data = df,
warmup = 1000, iter = 5000,
cores = parallel::detectCores()
)Distribution postérieure
post <- as_draws_df(x = mod5) # extracts posterior samples
R <- rethinking::rlkjcorr(n = 16000, K = 2, eta = 2) # samples from prior
data.frame(prior = R[, 1, 2], posterior = post$cor_cafe__Intercept__afternoon) %>%
gather(type, value, prior:posterior) %>%
ggplot(aes(x = value, color = type, fill = type) ) +
geom_histogram(position = "identity", alpha = 0.2) +
labs(x = expression(rho), y = "Nombre d'échantillons")
Shrinkage en deux dimensions
Comparaison de modèles
On compare le premier modèle (complete pooling model), le troisième modèle (partial pooling model), et le dernier modèle.
# comparaison des WAIC de chaque modèle
mod5 <- add_criterion(mod5, "waic")
w <- loo_compare(mod1, mod2, mod3, mod5, criterion = "waic")
print(w, simplify = FALSE) elpd_diff se_diff elpd_waic se_elpd_waic p_waic se_p_waic waic se_waic
mod5 0.0 0.0 -155.1 10.0 26.4 2.6 310.2 20.1
mod3 -98.9 8.3 -253.9 8.3 18.3 1.5 507.9 16.6
mod2 -99.8 8.3 -254.8 8.4 19.8 1.6 509.7 16.8
mod1 -156.0 13.7 -311.1 10.5 2.0 0.3 622.2 21.1 mod1 mod2 mod3 mod5
1.741945e-68 4.751919e-44 1.169499e-43 1.000000e+00 Comparaison de modèles
L’estimation du temps d’attente moyen est plus incertaine lorsqu’on prend en compte de nouvelles sources d’erreur. Cependant, l’erreur du modèle (i.e., ce qui n’est pas expliqué), la variation résiduelle \(\sigma\), diminue…
Estimate Est.Error Q2.5 Q97.5
b_Intercept 3.119875 0.07985599 2.960233 3.274415
sigma 1.143102 0.05749243 1.035762 1.263518 Estimate Est.Error Q2.5 Q97.5
b_Intercept 3.1192355 0.20255593 2.7183344 3.5238007
sigma 0.8222354 0.04359112 0.7414241 0.9121809 Estimate Est.Error Q2.5 Q97.5
b_Intercept 3.7346955 0.21471048 3.3059343 4.154400
b_afternoon -1.2334931 0.08701202 -1.4048471 -1.060672
sigma 0.4898649 0.02764033 0.4389085 0.547333Conclusions
Les modèles multi-niveaux (ou “modèles mixtes”) sont des extensions naturelles des modèles de régression classiques, où les paramètres de ces derniers se voient eux-même attribués des “modèles”, gouvernés par des hyper-paramètres.
Cette extension permet de faire des prédictions plus précises en prenant en compte la variabilité liée aux groupes ou structures (clusters) présent(e)s dans les données. Autrement dit, en modélisant les populations d’où sont tirés les effets aléatoires (e.g., la population de participants ou de stimuli).
Un modèle de régression classique est équivalent à un modèle multi-niveaux où la variabilité des effets aléatoires serait fixée à \(0\).
La cadre bayésien permet une interprétation naturelle des distributions desquelles proviennent les effets aléatoires (varying effects). En effet, ces distributions peuvent être interprétées comme des distributions a priori, dont les paramètres sont estimés à partir des données.
Travaux pratiques - sleepstudy
Reaction Days Subject
1 249.5600 0 308
2 258.7047 1 308
3 250.8006 2 308
4 321.4398 3 308
5 356.8519 4 308
6 414.6901 5 308
7 382.2038 6 308
8 290.1486 7 308
9 430.5853 8 308
10 466.3535 9 308
11 222.7339 0 309
12 205.2658 1 309
13 202.9778 2 309
14 204.7070 3 309
15 207.7161 4 309
16 215.9618 5 309
17 213.6303 6 309
18 217.7272 7 309
19 224.2957 8 309
20 237.3142 9 309Travaux pratiques - sleepstudy
Travaux pratiques - sleepstudy
À vous de construire les modèles mathématiques et les modèles brms correspondant aux modèles suivants :
- Modèle avec seulement l’effet fixe de
Days - Modèle avec l’effet fixe de
Days+ un effet aléatoire deSubject(varying intercept) - Modèle avec l’effet fixe de
Days+ un effet aléatoire deSubject(varying intercept) + un effet aléatoire deDays(varying slope)
Comparez ensuite ces modèles en utilisant les outils discutés aux cours précédents (e.g., WAIC) et concluez.
Proposition de solution
fmod0 <- lm(Reaction ~ Days, sleepstudy)
fmod1 <- lmer(Reaction ~ Days + (1|Subject), sleepstudy)
fmod2 <- lmer(Reaction ~ Days + (1 + Days|Subject), sleepstudy)
anova(fmod1, fmod2)Data: sleepstudy
Models:
fmod1: Reaction ~ Days + (1 | Subject)
fmod2: Reaction ~ Days + (1 + Days | Subject)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
fmod1 4 1802.1 1814.8 -897.04 1794.1
fmod2 6 1763.9 1783.1 -875.97 1751.9 42.139 2 7.072e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Proposition de solution
Proposition de solution
mod7 <- brm(
Reaction ~ 1 + Days + (1 | Subject),
prior = c(
set_prior("normal(200, 100)", class = "Intercept"),
set_prior("normal(0, 10)", class = "b"),
set_prior("cauchy(0, 10)", class = "sigma"),
set_prior("cauchy(0, 10)", class = "sd")
),
data = sleepstudy,
warmup = 1000, iter = 5000,
cores = parallel::detectCores()
)Proposition de solution
mod8 <- brm(
Reaction ~ 1 + Days + (1 + Days | Subject),
prior = c(
set_prior("normal(200, 100)", class = "Intercept"),
set_prior("normal(0, 10)", class = "b"),
set_prior("cauchy(0, 10)", class = "sigma"),
set_prior("cauchy(0, 10)", class = "sd")
),
data = sleepstudy,
warmup = 1000, iter = 5000,
cores = parallel::detectCores()
)Proposition de solution
# calcul du WAIC et ajout du WAIC à chaque modèle
mod6 <- add_criterion(mod6, "waic")
mod7 <- add_criterion(mod7, "waic")
mod8 <- add_criterion(mod8, "waic")
# comparaison des WAIC de chaque modèle
w <- loo_compare(mod6, mod7, mod8, criterion = "waic")
print(w, simplify = FALSE) elpd_diff se_diff elpd_waic se_elpd_waic p_waic se_p_waic waic se_waic
mod8 0.0 0.0 -860.5 22.3 32.9 8.3 1721.0 44.7
mod7 -24.3 11.5 -884.8 14.5 19.3 3.4 1769.6 28.9
mod6 -92.8 20.9 -953.3 10.6 3.2 0.5 1906.6 21.1 mod6 mod7 mod8
4.808945e-41 2.798348e-11 1.000000e+00 Références
Ladislas Nalborczyk - IMSB2022
Efron, B., & Morris, C. (1977). Stein’s paradox in statistics. Scientific American, 236(5), 119–127. https://doi.org/10.1038/scientificamerican0577-119
Gelman, A. (2005). Analysis of variance? Why it is more important than ever. The Annals of Statistics, 33(1), 1–53. https://doi.org/10.1214/009053604000001048
Gelman, A., & Hill, J. (2006). Data analysis using regression and multilevel/hierarchical models. https://doi.org/10.1017/cbo9780511790942
Lewandowski, D., Kurowicka, D., & Joe, H. (2009). Generating random correlation matrices based on vines and extended onion method. Journal of Multivariate Analysis, 100(9), 1989–2001. https://doi.org/10.1016/j.jmva.2009.04.008